Contrasts between coarsening and relaxational dynamics of surfaces 
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We discuss static and dynamic fluctuations of domain walls separating areas of constant but 
different slopes in steady-state configurations of crystalline surfaces both by an analytic treatment 
of the appropriate Langevin equation and by numerical simulations. In contrast to other situations 
that describe the dynamics in Ising-like systems such as models A and B, we find that the dynamic 
exponent 2 = 2 that governs the domain wall relaxation function is not equal to the inverse of the 
exponent n « 1/4 that describes the coarsening process that leads to the steady state. 

PACS numbers: 64.60. Cn, 05.70.Ln, 64.60.My, 68.35. Ja 



The description of a number of surface relaxation pro- 
cesses in which the primary relaxation mechanism is dif- 
fusion of particles along the interface begins, in a contin- 
uum picture, with a Langevin equation of the form 

a,Mi-,i)+V-j[Mr,i)]=77(r,t) (1) 

where r = (xi, X2) denotes the co-ordinates in a reference 
plane, h is the height of the surface above it, j represents 
a surface diffusion current, and rj accounts for thermal 
noise. The conservation law in (nl) reflects the assumed 
absence of voids, overhangs and evaporation of particles. 
For example, this equation was applied to the case of 
spinodal decomposition of an unstable crystal facet [|l],pl . 
In that case, V • j is simply the functional derivative of 
the surface free energy with respect to the function h. 

Our interest here is another case: the growth of films 
by molecular beam epitaxy (MBE). It was realized some 
time ago [^jji) that in MBE the breaking of detailed bal- 
ance by the deposition process and the presence of step- 
edge barriers could lead to a nonequilibrium diffusion cur- 
rent that destabilizes a singular surface and leads to the 
formation of mounds or pyramids ^j^ ■ This behavior is 
observed in simulations of discrete microscopic models as 
well as numerical integration of the following version of 
Eq. (I): 



dth{r, t) = -DAAh{r, t) - vAh{r, t) 



(2) 



where di, = djdxy, v = 1,2, and A = 9^ + 9|. The 
first term on the right arises from the equilibrium diffu- 
sion current [pj, the next two represent the nonequilib- 
rium diffusion current which with w, u > destabilize 
the singular surface {dyh = 0) and lead to a steady state 
composed of facets of slope dyh = ±^v/u. In partic- 
ular, during the late stages, both methods reveal that 
these three-dimensional features coarsen as function of 
time with the characteristic length scale R evolving as i" 
with 0.16 ^ n ^ 0.25. At present, this property is not 



well understood. 

In this letter, we focus on another aspect of the dy- 
namics of this model, namely, relaxation into a steady 
state consisting of two domains. We also investigate an 
intimately related subject: the fluctuations of the inter- 
face separating these domains. Here, a domain refers 
to an area of constant slope mentioned above. For the 
analytic part, we will rely on Eq. (y), setting, for sim- 
plicity, u — V — 2D = 1. To compare with these results 
we will present data from both simulations and numer- 
ical integration of (y). Our main conclusion is that z, 
the exponent associated with relaxational dynamics, is 
2. Specifically, we find that the relaxation time of a do- 
main wall fluctuation of wavevector k scales as Tk ~ k~'^ . 
Thus, we have a system which violates z — 1/n, a rela- 
tion considered to be generic, coming from the dynamic 
scaling assumption with only a single time scale. Two 
well-studied systems obeying this relation are models A 
and B ||l^. In the former, the AUen-Cahn equation [l^ ] 
used to describe coarsening dynamics is even identical 
to the equation which describes the slowest relaxational 
mode p4| , giving z — 1/n = 2. Similarly, the dynamics 
of model B is known to display a coarsening exponent 
n = i [ pl| and a relaxation exponent 2 = 3 |12|. In con- 
trast to this simple situation, the behavior of our system 
is considerably richer and more intriguing. Not only do 
we find z — 2 ^ 1/n, we also observe the steady-state 
structure factor to scale as S{k) ~ fc~^. These results in- 
dicate that the coarsening process is controlled by defects 
that are not present in the stationary state. 

It is tempting to compare our system with model B, 
both having similar non-linearities in the fields and con- 
servation laws. However, there are serious differences, 
which are most apparent if we rewrite \^ in terms of a 
two-component "order parameter", m^ = d^h: 

dtmy(r,t) = -iAAm^(r,t) - Am,.(r, i) 

+ ^a.94m3(r,i)]+9,7y(r,t), (3) 



with the constraint d^j^m^ = dj^rrifj,. By contrast, for a 
model B with two decoupled order parameters and a sim- 
ple free energy density !F = J2i, K^'^i/)^ ~ ^""^1 + j''^t^ 
the third term would read Am^, and no constraint would 
be imposed on ttt,^. 

We begin by discussing a stationary state with two do- 
main walls. Such a state can be imposed on a system with 
substrate dimensions Lx, Ly by screw boundary condi- 
tions in x: h{x + Lx,y) = h{x, y) + Lx and periodic ones 
in y: h{x, y + Ly) — h{x, y). Thus the average slope in x 
can be unity, which is one of the preferred ground states. 
In the y-direction, there is a single ridge and a single val- 
ley and, except for two regions of characteristic size i^ « 
^2D/v = 1, dyh{x,y) = ±1. We denote the steady-state 
height function by hc{x,y) = x + $c(y) and its deriva- 
tive with respect to y by ipdy) = dyhc{x,y) = ^'^y). 
Though <I>c can be found exactly, its details will turn out 
to be unimportant. To provide the reader some famil- 
iar ground, we quote an example in the limit ^/Ly — > 0: 
ifciy) -^ tanh{(j/ - l)/0 - tanh((y - 3l)/() - 1, with 
I = Ly/A, where the valley and ridge are located at y — I 
and y = M, respectively. 

Assuming now that h{r,t) = hc{r) + x(r, and lin- 
earizing (g), we find that x satisfies the Langevin equa- 
tion 

dtx{r,t) = -Cxir,t)+7jir,t) (4) 

with the Hcrmitian operator £ — ^AA—dx[V{dxhc)dx] — 
dy[V{dyhc)dy] aud V{C) = 3C - 1- The fluctuations of 
the surface can thus be expanded 



X 



ir,t)=J2Mt)viir) , 



(5) 



in terms of the eigenfunctions ve which solve the eigen- 
value problem 



Cvi{r) = XiVi{r) = X„{q)-)p„{q,y)e *«"= 



(6) 



Note that, for the last equality, we have exploited the 
translational invariance in x. 

There are several branches of the eigenvalue spectrum. 
Two of these, the height fluctuation modes and the 'Gold- 
stone' branch, terminate at q = with eigenvalue zero. 
The eigenfunctions of these zero-frequency modes are 
simply -ipni^iy) = 1 and 7/10(0,2/) = Vc{y)- The result 
^h{<1 = 0) = follows from translational invariance in 
h, while \Q[q = Q) — Q is due to Goldstone's theorem. 
The higher modes of the H-branch are simply obtained: 
ipniqjy) — exp{iqx) with Xniq) = 2(7^ + g'*/2. The non- 
trivial G- modes, on the other hand cannot be determined 
exactly. Focusing on small q, we write C — Ca+V{q) with 
V{q) = q^{2 — dy + q^) as a perturbation. Perturbation 
theory or the variational method 



X< dy ipCip / I dyip 



(7) 



with a trial function Lp 
0{q^lLy,q% 



Lpc both yield Xciq) = 2q^ 



Finally, a third set of low- lying modes can be identified 
— the breathing modes in which the ridge and valley 
move in opposite directions. An appropriate choice of the 
variational eigenfunction in (^ for the breathing mode is 
ipB{0,y) — tanh(j/—^)-|-tanh(j/—3Z)—?//Z-|-2 where the last 
two terms guarantee that the integral of this function is 
zero and that periodic boundary conditions are satisfied. 
A simple calculation using ip'siy) — IV'g(z^)I ~ 1/' then 
produces the bound Xsiq = 0) < 96/ Ly. It is also easy to 
see that higher modes have the property As (q) = c/Ly + 
2q^ +0{q'^/Ly, g"*) with some numerical constant c < 96. 
Thus, the breathing mode has a gap ~ l/-^y- However, 
for quadratic systems {Lx = Ly) this gap is of the same 
order as q^ and therefore Xb is of the same order as Xq 
for nonzero q. Thus, all eigenvalues that we have found 
so far scale in lowest order as A„ (g) ~ q'^ . For an exactly 
solvable version of this model ||l5[ it is found that this 
statement is indeed correct for all eigenvalues. Perhaps 
X ^ q^ should not be a surprise, considering the presence 
of at least two gradients on the RHS of (g) or (||). This is 
also the case, for a homogeneous state, in model B. The 
contrast with model B becomes dramatic only when we 
consider inhomogeneous states, which admit Goldstone 
modes. In model B, the fact that C is not Hermitian 
leads to A ~ q^ [|2| while here A ~ g^ . 

We now can calculate two-point correlation functions 
such as 5(k,r) = / d'^r'd'^r e*'"- (x(r', t)x{v' + r, i + r)) . 
Expanding the noise t] in terms of the eigenfunctions. 



we write 
Eqs. i) 
find 



2 T/f: 

J: 



■n 

T,t) ^ Lx^'^Y.nqBn{q,t)My)e-"'^- Using 
I) and (?7(r,t)?7(r',i)) = 2e5{v-v')5{t~t'), we 



S{Kt) 



E 



A«(g) 



V^4p)v5^4-p)e"^"^'^' 



(8) 



where k = (g,p) and ipnip) is the Fourier transform of 
the (normalized) ipniy)- This is a general result, show- 
ing that the term associated with the lowest eigenvalue 
will dominate the sum as t — > 00. However, in our case, 
there is no gap, i.e., all A's are 0{1/L'^), so that there 
is no valid reason to keep only one or two terms. Be- 
cause of the dependence on the Fourier transforms of the 
eigenfunctions, S{q^p,T) will not display simple scaling, 
apart from the special case p = 0. However, since all 
eigenvalues are proportional to fc^ the time-time corre- 
lation function <I>(k, r) = S (k, t) / S (k, 0) depends only 
on fc^T and thus the dynamic exponent that governs the 
fluctuations in the steady state is ^ = 2. 

We are also interested in the interface structure factor 
Si{q,T) as it may provide more insight into the coarsen- 
ing behavior. The position of the interfaces is determined 
by the condition dyh{x,y,t) — 0. This equation has two 
solutions yi{x,t) and y2{x,t) that correspond to the po- 
sitions of the top and the valley of the height profile. The 
interface structure factor can be calculated using either 
of the two functions, Si{q,t) — {yi{q,t + T)yi{—q,t)) = 
{y2{q,t + T)y2{—q,t)). In linear approximation we find 
that yi{x,t) = -dyx{v,t)/dlhc{y)\y=u and therefore 



Si{q,T) 



wm] 



2i9yX\y=iiQ,t + T)dyx\y^i{~q,t)) 



I^E[^.^"(0]^^' 



-A„(9)r 



(9) 



Because of the symmetries of the steady-state configura- 
tions all eigenfunctions are either symmetric or antisym- 
metric about y ~ and y — I. Of these four classes, only 
eigenfunctions that have a nonzero derivative a,t y — 1,31 
contribute to Sj{q,t) IM). The most important modes 
with this property are the Goldstone mode V'g(?/) ^-nd 
the breathing mode ipsiy)- Since dyijJG\y=h dyi^sly^i ^ 

— 1/2 

Ly we therefore find 



Si{q,T) = — 

Ijy 



^G -\G(q)i 



Xaiq) 



^B ^-\b{q)7 



Ai3(g) 



(10) 



where cq and cb are numeric constants. Hence the in- 
terfaces fluctuations are governed by a dynamic exponent 
z = 2 as well. Furthermore, Eq. ( [lO| ) shows that the in- 
terface structure factor scales as Si{q) ~ l/{q^Ly). A 
major consequence is that the interface width, u'(i) = 
[L~-^ ^ dxy\{x)Y/'^ ^ [Lx/Ly]^^^, does not diverge in the 
thermodynamic limit (provided L^ ~ Ly). Hence the 
pyramids that form in the coarsening dynamics are sep- 
arated by flat interfaces in striking contrast to the inter- 
faces in models A or B, in which the widths diverge as 

The coarsening dynamics of Eqs. (nh and (H) are still 
not well understood. Recently, Rost and Krug |16| pre- 
sented a scaling theory of coarsening for the rotationally 
invariant current jiso = ni(77iQ — m^) that neglects the 
crystalline symmetries of the growing film. It is ques- 
tionable whether such scaling theories are correct. The 
growth law for the average domain size depends usually 
on two lengthscales: the average domain size R{t) itself 
and the width of the domain walls ^. This was made clear 
by Bray and Rutenberg jlj] . For the Cahn-Hilliard equa- 
tion dtm = —A [Am + m{l ~ m^)] they find that the do- 
mains grow as R{t) ^ l/[^i?^(i)] leading to a coarsening 
exponent n = 1/3. Scaling and mean-field-like theories 
||l9| for the Cahn-Hilliard equation yield n = 1/4. This 
result is obtained, if the assumption is made that there is 
only one relevant lengthscale, the average domain size, so 
that ^ is replaced by R{t) in the expression above. Since 
the dimensionality of the operators in Eq. (@) is the same 
as those of the Cahn-Hilliard equation, scaling theories 
and dimensional analysis for equation IM) with a current 
j 7^ m/(|m|) are not very convincing. For the case of 
an isotropic current as studied in Rcf. n& the situation 
seems to be not so bad: At least in the analogous case of 
the Cahn-Hilliard equation only logarithmic corrections 
to the growth law are found p7[ if compared with scaling 
theories. 

Rost and Krug jl6| assume that the average slope 
approaches the limiting value ttiq in a power law man- 



ner, (m) = mo — ci~" , and they find that the expo- 
nent a' is linked to the coarsening exponent through 
a' = 2n. For model B the corresponding relation is 
different; the order parameter approaches its limiting 
value as toq — m{t) ~ t^^^^. Nevertheless, our calcu- 
lations do lend some support for the relation a' — 2n 
even in the case of an anisotropic surface current. A 
central result of our study is that all modes are gap- 
less, A(k) ~ k^. Consequently, slope fluctuations decay 
like (5m(t))2 ~ J d^r[d^xi^,t)]'^ - t'^/"^ and therefore 
a' = 1/2. Hence the relation a' = 2n yields n = 1/4 
in agreement with numerical solutions |B|,f8| of Eq. (0). 
Therefore, the differences between model B and surface 
dynamics that became apparent in the calculations pre- 
sented above, in particular the fact that n j^ z^^, may 
also be responsible for the difference in the growth law. 
We now discuss numerical calculations that verify the 
properties derived above. We have integrated Eq. (g) 
and also carried out Monte Carlo simulations for a dis- 
crete growth model ^^. Both display the same kind of 
instability and coarsening behavior. In the latter, parti- 
cles are randomly deposited on an initially flat substrate 
and the growing cluster relaxes by surface diffusion. An 
instability is produced by step-edge barriers in one of the 
two hopping directions on the square lattice. The steady 
state for a finite system is then a single mound with its 
ridge and valley perpendicular to this axis. In Fig. [^ we 
show the static structure factor Si{q) for the ridge for 



three different substrate sizes with L^ 



L. The 



collapse of the data, when Si{q) is multiplied by L is quite 
impressive. Further, the data seem to be converging to 
the (7~^-dependence predicted above. 

In Fig. H we show the relaxation functions ^{q,p,t) = 
Siq,p,t)/Siq,p,0) and $7(5, i) = Siiq,t)/Siiq,0) ob- 
tained by integration of Eq. (g) . The data collapse per- 
fectly to single curves when plotted as a function of 
q^t, p^t respectively, in complete agreement with (M) and 
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FIG. 1. The static structure factor for interface fluctua- 
tions in a model of MBE that is unstable toward mound for- 
mation Ml for substrate sizes L^ — Ly = 32, 64, 96 obtained 
from Monte Carlo simulations. Parameters are as given in 
pO[ . Note that the data collapse only if S{q) is multiplied by 
the substrate size L. 
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FIG. 2. Relaxation functions ${q,p,t) and $/(g, i) ob- 
tained from a numerical integration of Eq. (13) as a func- 
tion of the scaled arguments q^t and p^t, g„ = 2Trn/Lx and 



p„ = 2nn/Ly, L^ = Ly 



192. 



(p^). The only exceptions are the modes S{0,pn,t) with 
Pn — 2T:n/Ly and n an odd integer. For these modes, 
the difference of two large quantities is encountered and 
we have not succeeded in obtaining the necessary statis- 
tical and numerical accuracy to reach a definitive conclu- 
sion. More details of our study, including tests for the Ly 
dependence using rectangular lattices, will be presented 
elsewhere pSJ . 

In summary, we have demonstrated both by analytic 
calculations and by simulations that this class of nonequi- 
librium processes behaves quite differently from models A 
and B. In particular, the exponent z that controls the re- 
laxation of interface fluctuations in the steady state is not 
related to the coarsening exponent n by n = 1/z. While 
this is interesting in its own right, it also indicates that a 
theory of coarsening for these processes must be based on 
the interactions of more complicated defects than those 
allowed in the steady state. However, if the exponent re- 
lation a' — 2n ||lg is used, our result z = 2 = 1/a' yields 
a coarsening exponent n = 1/4. 

The results derived in this letter are also relevant for 
the interpretation of experimental results obtained, e.g., 
in scattering experiments of films grown using MBE. If 
the surface evolution can be described by pyramid for- 
mation and coarsening, one expects that the scattering 
amplitude has a "ring" structure ||21[] , i.e., a maximum 
at finite wavevectors kmaxlO ~ l/^(^) ~ *""• However, 
once fcmax becomes smaller than the minimum resolvable 
wavelength, the observed scattering is due to fluctuations 
of a single pyramid. In that case, the dynamic exponent 
z = 2 as explained above. Moreover, the measured rough- 
ness exponent a then also characterizes fluctuations su- 
perimposed on the pyramid structure and has the value 
(2 = 0. Experimentally, it may be very difficult to detect 
this logarithmic roughness. 
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